import cv2
from matplotlib import pyplot as plt
import numpy as np
import os
import scipy.signal as sgn
# Ivo's implemetation of Ellipse class
class Ellipse:
C : np.array #homogeneous conic matrix x'*C*x=0
def __init__( self, C : np.array ):
self.C = C
@classmethod
def estimate_ellipsoid( cls, points : np.array ):
"""
C = estimate_ellipsoid_2D ( points )
----------------------------------
estimates a homogeneous ellipsoid C
x'*C*x = 0
from a number of point samples
(in a least squares sense)
format of points
[ x^1_1 . . . . . x^N_1 ]
[ . . ]
points = [ . . ]
[ . . ]
[ x^1_n . . . . . x^N_n ]
[ 1 . . . . . 1 ]
where each point is element R^n
and there are N points, the points
have to be in homogeneous coordinates.
"""
n = points.shape[0]
print(n)
num_entries = int((n+1) * n/2)
if points.shape[1] < num_entries-1:
print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
# -------------- estimate ellipse ---------------
# assemble linear system
A = np.zeros( [points.shape[ 1 ], num_entries ] )
#for all points
for i in range( 0, points.shape[ 1 ] ):
#extract point
x = np.array([points[:, i]]).T
#print(x)
X = np.matmul( x, x.T )
#print(X)
X2 = X.T + X - np.diag( np.diag( X ) )
#print(X2)
row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]
#print(row)
A[ i, : ] = row
#print(A)
#solve for null space
U,D,V = np.linalg.svd( A )
ell = V[-1,:]
#create conic matrix
Cd = np.zeros( X2.shape )
Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell
C = Cd + Cd.T - np.diag( np.diag ( Cd ) )
#detC = np.linalg.det( C[0:-1,0:-1] )
#print(detC)
if Ellipse.is_hyperboloid( C ):
print('Warning: result is a hyperboloid!')
if Ellipse.is_paraboloid( C ):
print('Warning: result is a paraboloid!')
return C
@classmethod
def cov_mean_to_ellipsoid( cls, cov, mean, k ):
"""
C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from covariance matrix and mean value
of a gaussian distribution
k is a confidence factor determibed by the chi^2
distribution. The probability that x is inside
the ellipsoid is P_{chi^2}(k,2). Some values
for k are 2.41 (70%) and k=5.99 (95%).
cov : np.array - covariance matrix (non-homogeneous)
mean : np.array - mean value (vector)
k : float - confidence factor, see above
"""
#column vector, assume we get a 1D array
mean = np.array( [mean] ).T
SIGMAinv = np.linalg.inv( cov )
lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k
C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ),
np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
return C
@classmethod
def srt_to_ellipsoid( cls, S, R, T ):
"""
C = srt_to_ellipsoid_homog ( S, R, T )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from translation, rotation and scale parts as
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
"""
#function C = srt_to_ellipsoid_homog ( S, R, T )
invT = np.linalg.inv( T )
invR = R.T
invS = np.linalg.inv( S )
left = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
right = np.matmul( invS, np.matmul( invR, invT ) )
C = np.matmul( left, right )
return C
@classmethod
def decompose_ellipsoid( cls, C ):
"""
[S,R,T] = decompose_ellipsoid_homog ( C )
----------------------------------
decomposes a homogeneous ellipsoid C
x'*C*x = 0
into translation, rotation and scale parts
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
in some cases because the orientation of
the ellipsoid's major axis is not fixed.
The function generates a right handed system.
"""
#print(C)
C_ih = C[0:-1,0:-1]
#print(C_ih)
t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )
T = np.eye( C.shape[ 0 ] )
T[0:-1, -1 ] = t
U,D,V = np.linalg.svd( C_ih )
R = np.eye( C.shape[ 0 ] )
R[0:-1,0:-1] = V
#S = R' * T' * C * T * R;
S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )
S = S / S[-1,-1]
S = np.sqrt( np.diag(1/abs(np.diag(S))))
return S, R, T
@classmethod
def is_ellipsoid( cls, C ):
"""
bool = is_ellipsoid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is an ellipsoid
"""
#function val = conic_is_ellipsoid_homog ( C )
detC = np.linalg.det( C[0:-1,0:-1] )
return detC > 0
@classmethod
def is_hyperboloid( cls, C ):
"""
bool = is_hyperboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a hyperboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
return detC < 0
@classmethod
def is_paraboloid( cls, C ):
"""
bool = is_paraboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a paraboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
#maybe with EPS ?
return detC == 0
def get_center( self ):
S,R,T = self.decompose_ellipsoid( self.C )
return [ T[2,0], T[2,1] ]
def plot( self, ax, **kwargs ):
"""
plot_ellipsoid_2D ( C, style, N )
----------------------------------
plots a homogeneous ellipsoid C
x'*C*x = 0
style is a matlab plot style, e.g. 'b-' (default)
N is the number of polygon segments to
approximate the ellipse
"""
#function plot_ellipsoid_homog_2d ( C, style, N )
if not 'style' in locals():
style = 'b-'
if not 'N' in locals():
N = 100
angles = np.linspace( -np.pi, np.pi, N )
# generate a unit circle or a unit hyperbola
if ( Ellipse.is_ellipsoid( self.C ) ):
#uc = [ cos(x); sin(x); ones(size(x)) ];
x = np.cos( angles )
y = np.sin( angles )
w = np.ones( angles.shape )
uc = np.vstack([x,y,w])
else:
EPS = 0.1
x = np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
x2 = np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) ) #x2 = (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
uc = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
#decompose
S,R,T = Ellipse.decompose_ellipsoid( self.C )
#transform unit circle or parabola
transform = np.matmul( T, np.matmul( R, S ) )
ell = np.matmul( transform, uc )
if 'uc2' in locals():
ell2 = np.matmul( transform, uc2 )
ax.plot( ell2[0,:], ell2[1,:], style )
#plot conic
ax.plot( ell[0,:], ell[1,:], style )
#plot coordinate system
tx = T[0,2]
ty = T[1,2]
# print(tx, ty)
smin = S[0,0]
smaj = S[1,1]
ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )
ax.set_aspect('equal', adjustable='box')
return (tx, ty)
def fit_ellipse(path_dir_w_time, ifHull=False, date_now=""):
""" Find and fit an ellipse in an image.
it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
Parameters
----------
path_dir_w_time : str
_description_
date_now : str
_description_
Returns
-------
Ellipse
Ellipse in an image.
"""
I_temp = plt.imread(path_dir_w_time)[:,:,1].astype(float)/255
I = old_rplc_dead_px_med(I_temp)
ellipsoid_points=[]
if ifHull == False:
#every Nth scanline should be sampled
START_SCANLINE = 50 #there are some weird effects on the zero scanline that lead to detection of erroneous points
DELTA_SCANLINE = 100
#need to filter enough such that the big jump into and out of the disc has indeed the largest magnitude derivative
FILTERWIDTH=250
x=np.linspace(-4,4,FILTERWIDTH)
gauss=np.exp(-x**2)
plt.figure(figsize=[18,12])
for sl in range( START_SCANLINE, I.shape[0]-START_SCANLINE, DELTA_SCANLINE ):
#print(f'Checking Scanline {sl} ...' )
#line = self.I[sl]
line = I[sl, :]
filtered_line = sgn.convolve(line,gauss,mode='same')
#crude filter for line intersecting the ring structure
if np.amax(filtered_line) > 0.9:
#print("- scanline intersection ")
filtered_deriv = np.diff( filtered_line )
#crude extraction - this will probably need improvement
max_pos = np.argmax( filtered_deriv )
min_pos = np.argmin( filtered_deriv )
ellipsoid_points.append([max_pos,sl])
ellipsoid_points.append([min_pos,sl]) #add two points on perimeter
else:
#print("- no intersection ")
pass
print(ellipsoid_points)
else:
img = np.array(I*255, dtype=np.uint8)
cnt_img, final_cnt = contour_points(img)
if final_cnt is None:
return None, None
elif final_cnt.shape[0] < 5:
return None, None
else:
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
# print(ellipsoid_points)
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
print(ellipsoid_points)
#plt.clf()
#plt.imshow(self.I)
#plt.show()
for i in range( 0, ellipsoid_points.shape[1] ):
plt.plot(ellipsoid_points[0,i],ellipsoid_points[1,i],'rx')
E = Ellipse( Ellipse.estimate_ellipsoid( ellipsoid_points ) )
S,R,T=Ellipse.decompose_ellipsoid(E.C)
plt.plot(T[0,2],T[1,2],'gx')
plt.imshow(I)
center = E.plot(plt.gca())
plt.title(f'ellipsoide fit {1}')
# plt.savefig(f"{path_dir_w_time}/{date_now}_ledpos_{self.x_pos}_{self.y_pos}_{self.channel}_ellipsoide_fit.png")
return E, center
def old_rplc_dead_px_med(I_temp):
"""Replace dead pixels by median
Parameters
----------
I_temp : ndarray
image to be processed
Returns
-------
ndarray
clean image where dead pixels are replaced by median
"""
Imed = sgn.medfilt2d( I_temp, kernel_size=5 )
maxI = np.amax( Imed ) #normalize robustly
I_temp = I_temp / maxI
#replace dead pixels by median filtered values
I_temp[ I_temp > 1 ] = Imed[ I_temp > 1 ]
return I_temp
def contour_points(g_img):
kernel = np.ones((3, 3), np.uint8)
gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)
blur = cv2.GaussianBlur(gradient, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
gradient_img = np.zeros_like(g_img, dtype=np.uint8)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 2:
return None, None
cv2.drawContours(gradient_img, contours, -1, 255, 2)
kernel = np.ones((5, 5), np.uint8)
cv2.dilate(gradient_img, kernel, iterations = 1)
hull_img = np.zeros_like(g_img, dtype=np.uint8)
blur = cv2.GaussianBlur(g_img, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 2:
return None, None
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)
cv2.drawContours(hull_img, [hull], 0, 125, 2)
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
op[:10, :] = 0
op[-10:, :] = 0
kernel = np.ones((3, 3), np.uint8)
opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
x, y = np.where(opening > 1)
final_cnt = np.vstack([y, x]).T
# final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
imgs = np.hstack([gradient_img, hull_img, opening])
if final_cnt.shape[0] > 4:
return opening, final_cnt
else:
return None, None
dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_16_y_16_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
imgs, cnt = contour_points(g_img)
kernel = np.ones((25, 25),np.uint8)
cv2.dilate(imgs, kernel, iterations = 1)
fig, ax = plt.subplots(1, 1, figsize=(15, 15))
plt.imshow(imgs, cmap='gray')
plt.show()
ellipse = cv2.fitEllipse(cnt)
cv2.ellipse(g_img, ellipse, 255, 2)
fig, ax = plt.subplots(1, 1, figsize=(15, 15))
ax.imshow(g_img, cmap='gray')
plt.show()
print(ellipse)
((2971.643798828125, 913.9298095703125), (5347.255859375, 5739.8564453125), 169.41822814941406)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff"
# dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
# file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_16_y_16_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
g_img
array([[ 0, 0, 0, ..., 0, 0, 0],
[ 0, 0, 0, ..., 0, 0, 0],
[ 0, 0, 0, ..., 0, 0, 0],
...,
[ 887, 153, 1258, ..., 0, 0, 0],
[1531, 1343, 2342, ..., 0, 0, 0],
[1369, 1267, 1454, ..., 0, 0, 0]], dtype=uint16)
plt.imshow(g_img, cmap='gray')
<matplotlib.image.AxesImage at 0x208ad516f10>
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
plt.imshow(g_img_8, cmap='gray')
<matplotlib.image.AxesImage at 0x208b074b070>
kernel = np.ones((15, 15), np.uint8)
dilation = cv2.dilate(g_img_8, kernel, iterations = 1)
plt.imshow(dilation, cmap='gray')
plt.show()
kernel = np.ones((15, 15), np.uint8)
erosion = cv2.erode(g_img_8, kernel, iterations = 1)
plt.imshow(erosion, cmap='gray')
plt.show()
plt.imshow((np.fabs(-erosion.astype(np.int16) + dilation.astype(np.int16))).astype(np.uint8), cmap='gray')
plt.show()
kernel = np.ones((15, 15), np.uint8)
gradient = cv2.morphologyEx(g_img_8, cv2.MORPH_GRADIENT, kernel)
# g_img_2 = g_img.cop
plt.imshow(gradient, cmap='gray')
<matplotlib.image.AxesImage at 0x208b23ed430>
# kernel = np.ones((25, 25), np.uint8)
# opening = cv2.morphologyEx(gradient, cv2.MORPH_OPEN, kernel)
# plt.imshow(opening, cmap='gray')
blur = cv2.GaussianBlur(gradient+20, (5, 5), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
gradient_img = np.zeros_like(g_img, dtype=np.uint8)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 2:
print("Errors")
cv2.drawContours(gradient_img, contours, -1, 255, 2)
kernel = np.ones((35, 35),np.uint8)
gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)
plt.imshow(gradient_img, cmap='gray')
<matplotlib.image.AxesImage at 0x208d1f2adf0>
ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = gradient_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8+50
new_img = np.array((ch_img_1 + ch_img_2), dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
hull_img = np.zeros_like(g_img_8, dtype=np.uint8)
blur = cv2.GaussianBlur(g_img_8, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
plt.imshow(th3, cmap='gray')
<matplotlib.image.AxesImage at 0x208d0ee2070>
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)
if len(contours) < 2:
print("Errors")
cv2.drawContours(hull_img, [hull], 0, 125, 2)
plt.imshow(hull_img, cmap='gray')
<matplotlib.image.AxesImage at 0x208d0f92f10>
# hull_img = np.zeros_like(g_img_8, dtype=np.uint8)
# area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
# max_cntr_index = np.argmax(area_cntrs)
# cnt = contours[max_cntr_index]
# hull = cv2.convexHull(cnt)
# cv2.drawContours(hull_img, [hull], 0, 125, 2)
# plt.imshow(hull_img, cmap='gray')
ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = hull_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
op[:10, :] = 0
op[-10:, :] = 0
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(op, cmap='gray')
plt.show()
# op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
# op[:10, :] = 0
# op[-10:, :] = 0
ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = op
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8+50
new_img = np.array((ch_img_1 + ch_img_2), dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
# ellipse fit below not terrific but not because of contour above is v. good
x, y = np.where(op > 1)
op_cnt = np.vstack([y, x]).T
ellipse = cv2.fitEllipse(op_cnt)
cv2.ellipse(img, ellipse, (2**15, 0, 0), 3)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(cv2.convertScaleAbs(img, alpha=(255/65535)) + 125)
plt.show()
ellipse
((2975.01220703125, 1009.5517578125), (5253.75927734375, 5340.55322265625), 97.18701934814453)
min_axis_len = 5253.75927734375
maj_axis_len = 5340.55322265625
ratio = min_axis_len/maj_axis_len
ratio
0.9837481358777039
def contour_points_v2(g_img):
kernel = np.ones((15, 15), np.uint8)
gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)
blur = cv2.GaussianBlur(gradient + 20, (5, 5), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
gradient_img = np.zeros_like(g_img, dtype=np.uint8)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 1:
return None, None
cv2.drawContours(gradient_img, contours, -1, 255, 2)
kernel = np.ones((25, 25),np.uint8)
gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)
hull_img = np.zeros_like(g_img, dtype=np.uint8)
blur = cv2.GaussianBlur(g_img, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 1:
return None, None
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)
cv2.drawContours(hull_img, [hull], 0, 125, 2)
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
# op[:10, :] = 0
# op[-10:, :] = 0
# kernel = np.ones((3, 3), np.uint8)
# opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
x, y = np.where(op > 1)
final_cnt = np.vstack([y, x]).T
# final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
# imgs = np.hstack([gradient_img, hull_img, opening])
if final_cnt.shape[0] > 4:
return op, final_cnt
else:
return None, None
# now preprocessing fixed for both types of images
# Ivo's implemetation of Ellipse class
class Ellipse_Fixed_Ratio:
C : np.array #homogeneous conic matrix x'*C*x=0
def __init__( self, C : np.array ):
self.C = C
@classmethod
def estimate_ellipsoid( cls, points : np.array ):
"""
C = estimate_ellipsoid_2D ( points )
----------------------------------
estimates a homogeneous ellipsoid C
x'*C*x = 0
from a number of point samples
(in a least squares sense)
format of points
[ x^1_1 . . . . . x^N_1 ]
[ . . ]
points = [ . . ]
[ . . ]
[ x^1_n . . . . . x^N_n ]
[ 1 . . . . . 1 ]
where each point is element R^n
and there are N points, the points
have to be in homogeneous coordinates.
"""
n = points.shape[0]
# n is 3
num_entries = int((n+1) * n/2) - 2 # fixing ellipse maj/min axis length ratio
if points.shape[1] < num_entries-1:
print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
# -------------- estimate ellipse ---------------
# assemble linear system
A = np.zeros( [points.shape[ 1 ], num_entries ] )
#for all points
# print(points.shape[ 1 ])
for i in range( 0, points.shape[ 1 ] ):
#extract point
x = np.array([points[:, i]]).T
X = np.matmul( x, x.T )
X2 = X.T + X - np.diag( np.diag( X ) )
row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]
row[0] = row[0] + row[3]
row = np.delete(row, 3)
row = np.delete(row, 1)
A[ i, : ] = row
#solve for null space
U,D,V = np.linalg.svd( A )
# print(U, D, V)
ell = V[-1,:]
ell = np.insert(ell, 1, 0)
ell = np.insert(ell, 3, ell[0])
# ell[0] = 0.9837*ell[0]
print('Ell:', ell)
#create conic matrix
Cd = np.zeros( X2.shape )
Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell
C = Cd + Cd.T - np.diag( np.diag ( Cd ) )
print(f'C: {C}')
if Ellipse_Fixed_Ratio.is_hyperboloid( C ):
print('Warning: result is a hyperboloid!')
if Ellipse_Fixed_Ratio.is_paraboloid( C ):
print('Warning: result is a paraboloid!')
return C
@classmethod
def cov_mean_to_ellipsoid( cls, cov, mean, k ):
"""
C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from covariance matrix and mean value
of a gaussian distribution
k is a confidence factor determibed by the chi^2
distribution. The probability that x is inside
the ellipsoid is P_{chi^2}(k,2). Some values
for k are 2.41 (70%) and k=5.99 (95%).
cov : np.array - covariance matrix (non-homogeneous)
mean : np.array - mean value (vector)
k : float - confidence factor, see above
"""
#column vector, assume we get a 1D array
mean = np.array( [mean] ).T
SIGMAinv = np.linalg.inv( cov )
lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k
C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ),
np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
return C
@classmethod
def srt_to_ellipsoid( cls, S, R, T ):
"""
C = srt_to_ellipsoid_homog ( S, R, T )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from translation, rotation and scale parts as
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
"""
#function C = srt_to_ellipsoid_homog ( S, R, T )
invT = np.linalg.inv( T )
invR = R.T
invS = np.linalg.inv( S )
left = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
right = np.matmul( invS, np.matmul( invR, invT ) )
C = np.matmul( left, right )
return C
@classmethod
def decompose_ellipsoid( cls, C ):
"""
[S,R,T] = decompose_ellipsoid_homog ( C )
----------------------------------
decomposes a homogeneous ellipsoid C
x'*C*x = 0
into translation, rotation and scale parts
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
in some cases because the orientation of
the ellipsoid's major axis is not fixed.
The function generates a right handed system.
"""
#print(C)
C_ih = C[0:-1,0:-1]
#print(C_ih)
t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )
T = np.eye( C.shape[ 0 ] )
T[0:-1, -1 ] = t
U,D,V = np.linalg.svd( C_ih )
R = np.eye( C.shape[ 0 ] )
R[0:-1,0:-1] = V
#S = R' * T' * C * T * R;
S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )
S = S / S[-1,-1]
S = np.sqrt( np.diag(1/abs(np.diag(S))))
return S, R, T
@classmethod
def is_ellipsoid( cls, C ):
"""
bool = is_ellipsoid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is an ellipsoid
"""
#function val = conic_is_ellipsoid_homog ( C )
detC = np.linalg.det( C[0:-1,0:-1] )
return detC > 0
@classmethod
def is_hyperboloid( cls, C ):
"""
bool = is_hyperboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a hyperboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
return detC < 0
@classmethod
def is_paraboloid( cls, C ):
"""
bool = is_paraboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a paraboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
#maybe with EPS ?
return detC == 0
def get_center( self ):
S,R,T = self.decompose_ellipsoid( self.C )
return [ T[2,0], T[2,1] ]
def plot( self, ax, **kwargs ):
"""
plot_ellipsoid_2D ( C, style, N )
----------------------------------
plots a homogeneous ellipsoid C
x'*C*x = 0
style is a matlab plot style, e.g. 'b-' (default)
N is the number of polygon segments to
approximate the ellipse
"""
#function plot_ellipsoid_homog_2d ( C, style, N )
if not 'style' in locals():
style = 'b-'
if not 'N' in locals():
N = 100
angles = np.linspace( -np.pi, np.pi, N )
# generate a unit circle or a unit hyperbola
if ( Ellipse_Fixed_Ratio.is_ellipsoid( self.C ) ):
#uc = [ cos(x); sin(x); ones(size(x)) ];
x = np.cos( angles )
y = np.sin( angles )
w = np.ones( angles.shape )
uc = np.vstack([x,y,w])
else:
EPS = 0.1
x = np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
x2 = np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) ) #x2 = (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
uc = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
#decompose
S,R,T = Ellipse_Fixed_Ratio.decompose_ellipsoid( self.C )
#transform unit circle or parabola
transform = np.matmul( T, np.matmul( R, S ) )
ell = np.matmul( transform, uc )
if 'uc2' in locals():
ell2 = np.matmul( transform, uc2 )
ax.plot( ell2[0,:], ell2[1,:], style )
#plot conic
ax.plot( ell[0,:], ell[1,:], style )
#plot coordinate system
tx = T[0,2]
ty = T[1,2]
# print(tx, ty)
smin = S[0,0]
smaj = S[1,1]
ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )
ax.set_aspect('equal', adjustable='box')
return (tx, ty, smin, smaj)
def fit_ellipse_2(path_dir_w_time, ifHull=False, date_now=""):
""" Find and fit an ellipse in an image.
it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
Parameters
----------
path_dir_w_time : str
_description_
date_now : str
_description_
Returns
-------
Ellipse
Ellipse in an image.
"""
I_temp = plt.imread(path_dir_w_time)[:,:,1].astype(float)/255
I = old_rplc_dead_px_med(I_temp)
ellipsoid_points=[]
if ifHull == False:
#every Nth scanline should be sampled
START_SCANLINE = 50 #there are some weird effects on the zero scanline that lead to detection of erroneous points
DELTA_SCANLINE = 100
#need to filter enough such that the big jump into and out of the disc has indeed the largest magnitude derivative
FILTERWIDTH=250
x=np.linspace(-4,4,FILTERWIDTH)
gauss=np.exp(-x**2)
plt.figure(figsize=[18,12])
for sl in range( START_SCANLINE, I.shape[0]-START_SCANLINE, DELTA_SCANLINE ):
#print(f'Checking Scanline {sl} ...' )
#line = self.I[sl]
line = I[sl, :]
filtered_line = sgn.convolve(line,gauss,mode='same')
#crude filter for line intersecting the ring structure
if np.amax(filtered_line) > 0.9:
#print("- scanline intersection ")
filtered_deriv = np.diff( filtered_line )
#crude extraction - this will probably need improvement
max_pos = np.argmax( filtered_deriv )
min_pos = np.argmin( filtered_deriv )
ellipsoid_points.append([max_pos,sl])
ellipsoid_points.append([min_pos,sl]) #add two points on perimeter
else:
#print("- no intersection ")
pass
print(ellipsoid_points)
else:
img = np.array(I*255, dtype=np.uint8)
cnt_img, final_cnt = contour_points_v2(img)
if final_cnt is None:
return None, None
elif final_cnt.shape[0] < 5:
return None, None
else:
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
# print(ellipsoid_points)
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
# print(ellipsoid_points)
#plt.clf()
#plt.imshow(self.I)
#plt.show()
for i in range( 0, ellipsoid_points.shape[1] ):
plt.plot(ellipsoid_points[0,i],ellipsoid_points[1,i],'rx')
E = Ellipse_Fixed_Ratio( Ellipse_Fixed_Ratio.estimate_ellipsoid( ellipsoid_points ) )
S,R,T=Ellipse_Fixed_Ratio.decompose_ellipsoid(E.C)
print(f'S:{S}, \nR: {R}, \nT: {T}')
plt.plot(T[0,2],T[1,2],'gx')
plt.imshow(I)
center = E.plot(plt.gca())
plt.title(f'ellipsoide fit {1}')
# plt.savefig(f"{path_dir_w_time}/{date_now}_ledpos_{self.x_pos}_{self.y_pos}_{self.channel}_ellipsoide_fit.png")
return E, center
def old_rplc_dead_px_med(I_temp):
"""Replace dead pixels by median
Parameters
----------
I_temp : ndarray
image to be processed
Returns
-------
ndarray
clean image where dead pixels are replaced by median
"""
Imed = sgn.medfilt2d( I_temp, kernel_size=5 )
maxI = np.amax( Imed ) #normalize robustly
I_temp = I_temp / maxI
#replace dead pixels by median filtered values
I_temp[ I_temp > 1 ] = Imed[ I_temp > 1 ]
return I_temp
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"
# dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
# file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_15_y_15_r_0_g_1_b_0.tiff"
E, center = fit_ellipse(dir_path + file_path, ifHull=True)
plt.show()
print(center)
[[1.440e+02 1.450e+02 1.460e+02 1.440e+02 1.450e+02 1.460e+02 1.440e+02 1.450e+02 1.460e+02 4.557e+03 4.558e+03 4.559e+03 4.556e+03 4.557e+03 4.558e+03 4.559e+03 4.556e+03 4.557e+03 4.558e+03 4.559e+03 4.556e+03 4.557e+03 4.558e+03 4.480e+03 4.481e+03 4.482e+03 4.480e+03 4.481e+03 4.482e+03 4.480e+03 4.481e+03 4.482e+03 4.325e+03 4.326e+03 4.327e+03 4.325e+03 4.326e+03 4.327e+03 4.325e+03 4.326e+03 4.327e+03 3.509e+03 3.510e+03 3.511e+03 3.501e+03 3.502e+03 3.503e+03 3.504e+03 3.505e+03 3.509e+03 3.510e+03 3.511e+03 3.496e+03 3.497e+03 3.498e+03 3.499e+03 3.500e+03 3.501e+03 3.502e+03 3.503e+03 3.504e+03 3.505e+03 3.509e+03 3.510e+03 3.511e+03 3.492e+03 3.493e+03 3.494e+03 3.495e+03 3.496e+03 3.497e+03 3.498e+03 3.499e+03 3.500e+03 3.501e+03 3.502e+03 3.503e+03 3.504e+03 3.505e+03 3.486e+03 3.487e+03 3.488e+03 3.492e+03 3.493e+03 3.494e+03 3.495e+03 3.496e+03 3.497e+03 3.498e+03 3.499e+03 3.500e+03 3.486e+03 3.487e+03 3.488e+03 3.492e+03 3.493e+03 3.494e+03 3.495e+03 3.486e+03 3.487e+03 3.488e+03] [1.791e+03 1.791e+03 1.791e+03 1.792e+03 1.792e+03 1.792e+03 1.793e+03 1.793e+03 1.793e+03 3.192e+03 3.192e+03 3.192e+03 3.193e+03 3.193e+03 3.193e+03 3.193e+03 3.194e+03 3.194e+03 3.194e+03 3.194e+03 3.195e+03 3.195e+03 3.195e+03 3.265e+03 3.265e+03 3.265e+03 3.266e+03 3.266e+03 3.266e+03 3.267e+03 3.267e+03 3.267e+03 3.367e+03 3.367e+03 3.367e+03 3.368e+03 3.368e+03 3.368e+03 3.369e+03 3.369e+03 3.369e+03 3.730e+03 3.730e+03 3.730e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.736e+03 3.736e+03 3.736e+03] [1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00]] 3 Warning: result is a hyperboloid!
(5272.143308579376, 3112.6301639181233)
E, center = fit_ellipse_2(dir_path + file_path, ifHull=True)
plt.show()
print(center)
Ell: [ 5.23202987e-07 0.00000000e+00 -1.44310664e-03 5.14674778e-07 -6.28543710e-04 9.99998761e-01] C: [[ 5.23202987e-07 0.00000000e+00 -1.44310664e-03] [ 0.00000000e+00 5.14674778e-07 -6.28543710e-04] [-1.44310664e-03 -6.28543710e-04 9.99998761e-01]] S:[[2.67648633e+03 0.00000000e+00 0.00000000e+00] [0.00000000e+00 2.69857004e+03 0.00000000e+00] [0.00000000e+00 0.00000000e+00 1.00000000e+00]], R: [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]], T: [[1.00000000e+00 0.00000000e+00 2.75821559e+03] [0.00000000e+00 1.00000000e+00 1.22124444e+03] [0.00000000e+00 0.00000000e+00 1.00000000e+00]]
(2758.2155895125884, 1221.2444375066455, 2676.486329255531, 2698.5700358857857)
center[2]/center[3]
0.9918165152889924
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, 195, axis=1))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(20, 20))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_19_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -195), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -195), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -210), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_19_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -200), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (177, -180), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (187, -205), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img_2, _ = cv2.split(img)
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img
ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2
new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, -200), axis=(1, 0)))*2, dtype=np.uint8)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path)
_, g_img, _ = cv2.split(img)
g_img_2 = g_img.copy()
ellipse = ((2755+390+15, 1210-390),
(5322, 5200),
0)
# 17,17 -> not a good fit: more elliptical
# ((2755+390+15, 1210-390),
# (5322, 5200),
# 0)
# 13,13
# ((2755-390, 1210+390),
# (5322, 5322),
# 0)
# 15, 17
# ((2755, 1210-380),
# (5322, 5200),
# 0)
# 15, 13
# ((2755, 1210+400/390),
# (5320, 5320),
# 0)
# 13,15:
# ((2765-410, 1210-15),
# (5320, 5320),
# 0)
# 15,15:
# ((2765, 1210),
# (5320, 5320),
# 0)
# 17,15:
# ((2765+390, 1210-20),
# (5320, 5320),
# 0)
cv2.ellipse(g_img_2, ellipse, 255, 4)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(g_img_2, cmap='gray')
plt.show()
def contour_points_v2(g_img):
kernel = np.ones((15, 15), np.uint8)
gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)
blur = cv2.GaussianBlur(gradient + 20, (5, 5), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
gradient_img = np.zeros_like(g_img, dtype=np.uint8)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 1:
return None
cv2.drawContours(gradient_img, contours, -1, 255, 2)
kernel = np.ones((25, 25),np.uint8)
gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)
hull_img = np.zeros_like(g_img, dtype=np.uint8)
blur = cv2.GaussianBlur(g_img, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 1:
return None
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)
cv2.drawContours(hull_img, [hull], 0, 125, 2)
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255
# op[:10, :] = 0
# op[-10:, :] = 0
# kernel = np.ones((3, 3), np.uint8)
# opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
# final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
# imgs = np.hstack([gradient_img, hull_img, opening])
# if final_cnt.shape[0] > 4:
# return op
# else:
# return None
return op
# Ivo's implemetation of Ellipse class
class Ellipse:
C : np.array #homogeneous conic matrix x'*C*x=0
def __init__( self, C : np.array ):
self.C = C
@classmethod
def estimate_ellipsoid( cls, points : np.array ):
"""
C = estimate_ellipsoid_2D ( points )
----------------------------------
estimates a homogeneous ellipsoid C
x'*C*x = 0
from a number of point samples
(in a least squares sense)
format of points
[ x^1_1 . . . . . x^N_1 ]
[ . . ]
points = [ . . ]
[ . . ]
[ x^1_n . . . . . x^N_n ]
[ 1 . . . . . 1 ]
where each point is element R^n
and there are N points, the points
have to be in homogeneous coordinates.
"""
n = points.shape[0]
num_entries = int((n+1) * n/2)
if points.shape[1] < num_entries-1:
print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
# -------------- estimate ellipse ---------------
# assemble linear system
A = np.zeros( [points.shape[ 1 ], num_entries ] )
#for all points
for i in range( 0, points.shape[ 1 ] ):
#extract point
x = np.array([points[:, i]]).T
#print(x)
X = np.matmul( x, x.T )
#print(X)
X2 = X.T + X - np.diag( np.diag( X ) )
#print(X2)
row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]
#print(row)
A[ i, : ] = row
#print(A)
#solve for null space
U,D,V = np.linalg.svd( A )
ell = V[-1,:]
#create conic matrix
Cd = np.zeros( X2.shape )
Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell
C = Cd + Cd.T - np.diag( np.diag ( Cd ) )
#detC = np.linalg.det( C[0:-1,0:-1] )
#print(detC)
if Ellipse.is_hyperboloid( C ):
print('Warning: result is a hyperboloid!')
if Ellipse.is_paraboloid( C ):
print('Warning: result is a paraboloid!')
return C
@classmethod
def cov_mean_to_ellipsoid( cls, cov, mean, k ):
"""
C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from covariance matrix and mean value
of a gaussian distribution
k is a confidence factor determibed by the chi^2
distribution. The probability that x is inside
the ellipsoid is P_{chi^2}(k,2). Some values
for k are 2.41 (70%) and k=5.99 (95%).
cov : np.array - covariance matrix (non-homogeneous)
mean : np.array - mean value (vector)
k : float - confidence factor, see above
"""
#column vector, assume we get a 1D array
mean = np.array( [mean] ).T
SIGMAinv = np.linalg.inv( cov )
lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k
C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ),
np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
return C
@classmethod
def srt_to_ellipsoid( cls, S, R, T ):
"""
C = srt_to_ellipsoid_homog ( S, R, T )
----------------------------------
generates a homogeneous ellipsoid C
x'*C*x = 0
from translation, rotation and scale parts as
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
"""
#function C = srt_to_ellipsoid_homog ( S, R, T )
invT = np.linalg.inv( T )
invR = R.T
invS = np.linalg.inv( S )
left = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
right = np.matmul( invS, np.matmul( invR, invT ) )
C = np.matmul( left, right )
return C
@classmethod
def decompose_ellipsoid( cls, C ):
"""
[S,R,T] = decompose_ellipsoid_homog ( C )
----------------------------------
decomposes a homogeneous ellipsoid C
x'*C*x = 0
into translation, rotation and scale parts
x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
in some cases because the orientation of
the ellipsoid's major axis is not fixed.
The function generates a right handed system.
"""
#print(C)
C_ih = C[0:-1,0:-1]
#print(C_ih)
t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )
T = np.eye( C.shape[ 0 ] )
T[0:-1, -1 ] = t
U,D,V = np.linalg.svd( C_ih )
R = np.eye( C.shape[ 0 ] )
R[0:-1,0:-1] = V
#S = R' * T' * C * T * R;
S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )
S = S / S[-1,-1]
S = np.sqrt( np.diag(1/abs(np.diag(S))))
return S, R, T
@classmethod
def is_ellipsoid( cls, C ):
"""
bool = is_ellipsoid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is an ellipsoid
"""
#function val = conic_is_ellipsoid_homog ( C )
detC = np.linalg.det( C[0:-1,0:-1] )
return detC > 0
@classmethod
def is_hyperboloid( cls, C ):
"""
bool = is_hyperboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a hyperboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
return detC < 0
@classmethod
def is_paraboloid( cls, C ):
"""
bool = is_paraboloid_homog ( C )
----------------------------------
determines whether homogeneous conic C
x'*C*x = 0
is a paraboloid
"""
detC = np.linalg.det( C[0:-1,0:-1] )
#maybe with EPS ?
return detC == 0
def get_center( self ):
S,R,T = self.decompose_ellipsoid( self.C )
return [ T[2,0], T[2,1] ]
def plot( self, ax, **kwargs ):
"""
plot_ellipsoid_2D ( C, style, N )
----------------------------------
plots a homogeneous ellipsoid C
x'*C*x = 0
style is a matlab plot style, e.g. 'b-' (default)
N is the number of polygon segments to
approximate the ellipse
"""
#function plot_ellipsoid_homog_2d ( C, style, N )
if not 'style' in locals():
style = 'b-'
if not 'N' in locals():
N = 100
angles = np.linspace( -np.pi, np.pi, N )
# generate a unit circle or a unit hyperbola
if ( Ellipse.is_ellipsoid( self.C ) ):
#uc = [ cos(x); sin(x); ones(size(x)) ];
x = np.cos( angles )
y = np.sin( angles )
w = np.ones( angles.shape )
uc = np.vstack([x,y,w])
else:
EPS = 0.1
x = np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
x2 = np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) ) #x2 = (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
uc = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
#decompose
S,R,T = Ellipse.decompose_ellipsoid( self.C )
#transform unit circle or parabola
transform = np.matmul( T, np.matmul( R, S ) )
ell = np.matmul( transform, uc )
if 'uc2' in locals():
ell2 = np.matmul( transform, uc2 )
ax.plot( ell2[0,:], ell2[1,:], style )
#plot conic
ax.plot( ell[0,:], ell[1,:], style )
#plot coordinate system
tx = T[0,2]
ty = T[1,2]
# print(tx, ty)
smin = S[0,0]
smaj = S[1,1]
ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )
ax.set_aspect('equal', adjustable='box')
return (tx, ty)
def fit_ellipse_3(img_name, img_8, g_img_8, ifHull=False, date_now="", method="ivo", save_path=None):
""" Find and fit an ellipse in an image.
it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
Parameters
----------
path_dir_w_time : str
_description_
date_now : str
_description_
Returns
-------
Ellipse
Ellipse in an image.
"""
ellipsoid_points=[]
img = np.array(g_img_8, dtype=np.uint8)
op = contour_points_v2(img)
op[:20, :] = 0
op[-20:, :] = 0
op[:, :20] = 0
op[:, -20:] = 0
# plt.imshow(op)
# plt.show()
x, y = np.where(op > 1)
final_cnt = np.vstack([y, x]).T
if final_cnt.shape[0] > 10000:
final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)
if final_cnt is None:
return None, None, None
else:
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
# print(ellipsoid_points)
is_det_ellipse = True
if method == "ivo":
E = Ellipse( Ellipse.estimate_ellipsoid( ellipsoid_points ) )
if Ellipse.is_ellipsoid( E.C ):
S,R,T = Ellipse.decompose_ellipsoid(E.C)
tx = T[0,2]
ty = T[1,2]
smin = S[0,0]
smaj = S[1,1]
# f = E.C[-1, -1]
f = np.arctan2(R[1,0], R[0, 0])
# print(theta)
# print(E.C)
else:
is_det_ellipse = False
else:
(tx, ty), (smin, smaj), f = cv2.fitEllipse(ellipsoid_points[:2, :].T.astype(np.int32))
smin, smaj = smin/2, smaj/2
# print((tx, ty), (smin, smaj), f)
#
if is_det_ellipse:
for i in range( 0, ellipsoid_points.shape[1] ):
draw_cross(img_8, (ellipsoid_points[0,i],ellipsoid_points[1,i]), (0, 0, 255), length=20, thickness=5)
cv2.ellipse(img_8, (int(tx), int(ty)), (int(smin), int(smaj)), f, 0, 360, (255, 0, 0), 5)
draw_cross(img_8, (tx, ty), (0, 255, 0))
# plt.imshow(img_8)
# plt.show()
cv2.imshow('Resized Image', cv2.resize(img_8, (1080, 720), interpolation = cv2.INTER_LINEAR))
cv2.waitKey(0)
cv2.destroyAllWindows()
if save_path is not None:
cv2.imwrite(f'{save_path}/{img_name[:-5]}.png', img_8)
return (int(tx), int(ty)), (int(smin), int(smaj)), f
else:
return None, None, None
cv2.fitEllipse(ellipsoid_points[:2, :].T.astype(np.int32))
((219.3291473388672, 104.70303344726562), (513.478759765625, 532.2569580078125), 53.40561294555664)
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
center, axis, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='ivo')
[[ 7.23449895e-07 7.32307094e-09 -1.69863270e-03] [ 7.32307094e-09 7.46098008e-07 -4.49513750e-04] [-1.69863270e-03 -4.49513750e-04 -9.99998456e-01]]
def draw_cross(img, cord, color, length=50, thickness=10):
x, y = int(cord[0]), int(cord[1])
cv2.line(img, (x-length, y-length), (x+length, y+length), color, thickness)
cv2.line(img, (x+length, y-length), (x-length, y+length), color, thickness)
return img
# # grid based on thumbnail images
# fit_range = np.zeros([7, 7], dtype=np.uint8)
# center_cords_ivo = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_ivo = np.zeros([7, 7, 2], dtype=np.float64)
# for x in range(12, 19):
# for y in range(12, 19):
# dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
# file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
# print('x:', x, 'y:', y, file_path)
# img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
# _, g_img, _ = cv2.split(img)
# img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
# g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
# center, axis_size, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True)
# if center is not None:
# good_fit = input("Is it a good fit?")
# fit_range[x-12, y-12] = int(good_fit)
# center_cords_ivo[x-12, y-12] = center
# axis_lens_ivo[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff Is it a good fit?5 x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff Is it a good fit?3 x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff Is it a good fit?4 x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff Is it a good fit?2 x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff Is it a good fit?2 x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff Warning: result is a hyperboloid! x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff Is it a good fit?3 x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff Is it a good fit?3 x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff Warning: result is a hyperboloid! x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff Is it a good fit?2 x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff Is it a good fit?3 x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff Is it a good fit?3 x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff Warning: result is a paraboloid! x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff Is it a good fit?2 x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff Warning: result is a paraboloid! x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff Is it a good fit?3 x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff Is it a good fit?2 x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff Is it a good fit?3 x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff Is it a good fit?3 x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff Warning: result is a paraboloid! x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff Is it a good fit?3 x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff Is it a good fit?3 x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff Is it a good fit?5 x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff Warning: result is a paraboloid! x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff Is it a good fit?3 x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff Is it a good fit?2 x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff Is it a good fit?5
plt.figure(figsize=(15, 10))
center_leds = center_cords_ivo[3, 3]
y_unit_len = center_cords_ivo[3, 2][1] - center_cords_ivo[3, 3][1]
x_unit_len = center_cords_ivo[4, 3][0] - center_cords_ivo[3, 3][0]
center_leds, y_unit_len, x_unit_len
x_axis = [center_leds[0] - i*x_unit_len for i in range(-6, 6)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-5, 5)]
xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')
x_cords = [x if x != 0.0 else np.nan for x in center_cords_ivo.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_ivo.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range.reshape(-1)]
plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()
count = 0
for x in range(12, 19):
for y in range(12, 19):
plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
count += 1
plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
axis_lens_pro = axis_lens_ivo.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')
xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])
cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')
cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')
plt.show()
axis_lens_pro = axis_lens_ivo.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
# # grid based on thumbnail images
# fit_range_cv2 = np.zeros([7, 7], dtype=np.uint8)
# center_cords_cv2 = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_cv2 = np.zeros([7, 7, 2], dtype=np.float64)
# save_dir = "C:/Users/Kazim/Desktop/Overall/VignettingCircleDetection/310823_ellipse_cv2_fit_tiff_size"
# for x in range(12, 19):
# for y in range(12, 19):
# dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
# file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
# print('x:', x, 'y:', y, file_path)
# img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
# _, g_img, _ = cv2.split(img)
# img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
# g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
# center, axis_size, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='opencv', save_path=save_dir)
# if center is not None:
# good_fit = input("Is it a good fit?")
# fit_range_cv2[x-12, y-12] = int(good_fit)
# center_cords_cv2[x-12, y-12] = center
# axis_lens_cv2[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff Is it a good fit?5 x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff Is it a good fit?2 x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff Is it a good fit?2 x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff Is it a good fit?5 x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff Is it a good fit?2 x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff Is it a good fit?1 x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff Is it a good fit?5
plt.figure(figsize=(15, 10))
center_leds = center_cords_cv2[3, 3]
y_unit_len = center_cords_cv2[3, 2][1] - center_cords_cv2[3, 3][1]
x_unit_len = center_cords_cv2[4, 3][0] - center_cords_cv2[3, 3][0]
center_leds, y_unit_len, x_unit_len
x_axis = [center_leds[0] - i*x_unit_len for i in range(-4, 4)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-4, 4)]
xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')
x_cords = [x if x != 0.0 else np.nan for x in center_cords_cv2.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_cv2.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range_cv2.reshape(-1)]
plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()
count = 0
for x in range(12, 19):
for y in range(12, 19):
plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
count += 1
plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
axis_lens_pro = axis_lens_cv2.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')
xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])
cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')
cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')
plt.show()
axis_lens_pro = axis_lens_cv2.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
# fixing A, B, C, F
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
center, axis, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='ivo') # good fit
[[-3.24326055e-07 4.75306325e-10 9.59646019e-04] [ 4.75306325e-10 -3.27813554e-07 3.93645616e-04] [ 9.59646019e-04 3.93645616e-04 -9.99999462e-01]]
np.cos()
img = np.array(g_img_8, dtype=np.uint8)
op = contour_points_v2(img)
op[:20, :] = 0
op[-20:, :] = 0
op[:, :20] = 0
op[:, -20:] = 0
x, y = np.where(op > 1)
final_cnt = np.vstack([y, x]).T
if final_cnt.shape[0] > 10000:
final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points
array([[ 529., 530., 531., ..., 2178., 2179., 2180.],
[ 93., 93., 93., ..., 3732., 3732., 3732.]])
b = -3.24326055e-07*x**2 + -3.27813554e-07*y**2 + f
A = np.array([x, y]).T
A, -b
(array([[ 93, 529],
[ 93, 530],
[ 93, 531],
...,
[3732, 2178],
[3732, 2179],
[3732, 2180]], dtype=int64),
array([1.09454023, 1.09488739, 1.0952352 , ..., 7.07219959, 7.07362787,
7.07505681]))
np.linalg.lstsq(A, -b, rcond=None)
(array([0.00077151, 0.00194168]), array([0.65526845]), 2, array([338927.99799414, 121566.05420793]))
9.59646019e-04*2, 3.93645616e-04*2
(0.001919292038, 0.000787291232)
9.26610761e-04*2, -9.39963645e-04*2
(0.001853221522, -0.00187992729)
b = 4.75306325e-10*x*y + -9.99999462e-01
A = np.array([x**2, y**2, x, y]).T
A, -b
(array([[ 8649, 279841, 93, 529],
[ 8649, 280900, 93, 530],
[ 8649, 281961, 93, 531],
...,
[13927824, 4743684, 3732, 2178],
[13927824, 4748041, 3732, 2179],
[13927824, 4752400, 3732, 2180]], dtype=int64),
array([0.99997608, 0.99997603, 0.99997599, ..., 0.99613603, 0.99613426,
0.99613248]))
A, B, D, E = np.linalg.lstsq(A, -b, rcond=None)[0]
A, B, D/2, E/2
(-3.268969844369404e-07, -3.236289044146954e-07, 0.0003930310412070922, 0.0009580904040238742)
def fit_ellipse_4(img_name, img_8, g_img_8, save_path=None):
""" Find ellipse center position vector where A=-3.24326055e-07, B=0, C=-3.27813554e-07, F=-9.99999462e-01
Parameters
----------
Returns
-------
Ellipse
Ellipse in an image.
"""
ellipsoid_points=[]
img = np.array(g_img_8, dtype=np.uint8)
op = contour_points_v2(img)
op[:20, :] = 0
op[-20:, :] = 0
op[:, :20] = 0
op[:, -20:] = 0
x, y = np.where(op > 0)
final_cnt = np.vstack([y, x]).T
if final_cnt is None:
return None, None, None
else:
if final_cnt.shape[0] > 10000:
final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
is_det_ellipse = True
# C = -3.24326055e-07
B = 4.75306325e-10
# A = -3.27813554e-07
F = -9.99999462e-01*2
# b = A*x**2 + B*x*y + C*y**2 + F
b = B*x*y + F
# eqs = np.array([x**2, y**2, x, y]).T
# C, A, E, D = np.linalg.lstsq(eqs, -b, rcond=None)[0]/2
# matC = np.array([[A, B, D/2], [B, C, E/2], [D/2, E/2, F]])
eqs = np.array([x**2 + y**2, x, y]).T
A, E, D = np.linalg.lstsq(eqs, -b, rcond=None)[0]/2
matC = np.array([[A, B, D/2], [B, A, E/2], [D/2, E/2, F/2]])
E = Ellipse( matC )
# print(Ellipse.is_ellipsoid( E.C ))
print(matC)
# print(np.linalg.lstsq(eqs, -b, rcond=None))
if Ellipse.is_ellipsoid( E.C ):
S,R,T = Ellipse.decompose_ellipsoid(E.C)
tx = T[0,2]
ty = T[1,2]
smin = S[0,0]
smaj = S[1,1]
theta = np.arctan2(R[1,0], R[0, 0])
# print(theta)
# f = E.C[-1, -1]
# print(E.C)
else:
is_det_ellipse = False
if is_det_ellipse:
for i in range( 0, ellipsoid_points.shape[1] ):
draw_cross(img_8, (ellipsoid_points[0,i],ellipsoid_points[1,i]), (0, 0, 255), length=20, thickness=5)
cv2.ellipse(img_8, (int(tx), int(ty)), (int(smin), int(smaj)), theta, 0, 360, (255, 0, 0), 5)
draw_cross(img_8, (tx, ty), (0, 255, 0))
# plt.imshow(img_8)
# plt.show()
cv2.imshow('Resized Image', cv2.resize(img_8, (1080, 720), interpolation = cv2.INTER_LINEAR))
cv2.waitKey(0)
cv2.destroyAllWindows()
if save_path is not None:
cv2.imwrite(f'{save_path}/{img_name[:-5]}.png', img_8)
return (int(tx), int(ty)), (int(smin), int(smaj)), theta
else:
return None, None, None
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
center, axis_size, f = fit_ellipse_4(file_path, img_8, g_img_8, save_path=None)
[[-9.96262090e-07 4.75306325e-10 2.74601838e-03] [ 4.75306325e-10 -9.96262090e-07 7.62551767e-04] [ 2.74601838e-03 7.62551767e-04 -9.99999462e-01]]
# # grid based on thumbnail images
# fit_range_circle = np.zeros([7, 7], dtype=np.uint8)
# center_cords_circle = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_circle = np.zeros([7, 7, 2], dtype=np.float64)
# save_dir = "C:/Users/Kazim/Desktop/Overall/VignettingCircleDetection/310823_circle_fit_tiff_size"
# for x in range(12, 19):
# for y in range(12, 19):
# dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
# file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
# print('x:', x, 'y:', y, file_path)
# img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
# _, g_img, _ = cv2.split(img)
# img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
# g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
# center, axis_size, f = fit_ellipse_4(file_path, img_8, g_img_8, save_path=save_dir)
# if center is not None:
# good_fit = input("Is it a good fit?")
# fit_range_circle[x-12, y-12] = int(good_fit)
# center_cords_circle[x-12, y-12] = center
# axis_lens_circle[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff [[-1.51944827e-07 4.75306325e-10 4.28640374e-04] [ 4.75306325e-10 -1.51944827e-07 2.64667160e-04] [ 4.28640374e-04 2.64667160e-04 -9.99999462e-01]] Is it a good fit?5 x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff [[-8.05421194e-07 4.75306325e-10 1.84482464e-03] [ 4.75306325e-10 -8.05421194e-07 1.29403723e-03] [ 1.84482464e-03 1.29403723e-03 -9.99999462e-01]] Is it a good fit?3 x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff [[ 1.82876905e-06 4.75306325e-10 -3.95544562e-03] [ 4.75306325e-10 1.82876905e-06 -2.53086219e-03] [-3.95544562e-03 -2.53086219e-03 -9.99999462e-01]] Is it a good fit?1 x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff [[ 8.30978298e-07 4.75306325e-10 -1.78516376e-03] [ 4.75306325e-10 8.30978298e-07 -9.72636790e-04] [-1.78516376e-03 -9.72636790e-04 -9.99999462e-01]] Is it a good fit?2 x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff [[ 5.53879588e-07 4.75306325e-10 -1.18126017e-03] [ 4.75306325e-10 5.53879588e-07 -5.30531858e-04] [-1.18126017e-03 -5.30531858e-04 -9.99999462e-01]] Is it a good fit?2 x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff [[ 4.35008066e-07 4.75306325e-10 -9.22324007e-04] [ 4.75306325e-10 4.35008066e-07 -3.19552659e-04] [-9.22324007e-04 -3.19552659e-04 -9.99999462e-01]] Is it a good fit?2 x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff [[ 3.82984548e-07 4.75306325e-10 -8.12447687e-04] [ 4.75306325e-10 3.82984548e-07 -1.98597659e-04] [-8.12447687e-04 -1.98597659e-04 -9.99999462e-01]] Is it a good fit?2 x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff [[-5.58856538e-07 4.75306325e-10 1.33007273e-03] [ 4.75306325e-10 -5.58856538e-07 1.00255279e-03] [ 1.33007273e-03 1.00255279e-03 -9.99999462e-01]] Is it a good fit?1 x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff [[-9.23869595e-07 4.75306325e-10 2.18874479e-03] [ 4.75306325e-10 -9.23869595e-07 1.48275359e-03] [ 2.18874479e-03 1.48275359e-03 -9.99999462e-01]] Is it a good fit?1 x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff [[-2.27036215e-06 4.75306325e-10 5.37034775e-03] [ 4.75306325e-10 -2.27036215e-06 3.18468554e-03] [ 5.37034775e-03 3.18468554e-03 -9.99999462e-01]] Is it a good fit?1 x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff [[ 3.06726902e-06 4.75306325e-10 -7.19762759e-03] [ 4.75306325e-10 3.06726902e-06 -3.60634197e-03] [-7.19762759e-03 -3.60634197e-03 -9.99999462e-01]] Is it a good fit?2 x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff [[ 1.10581527e-06 4.75306325e-10 -2.58177670e-03] [ 4.75306325e-10 1.10581527e-06 -1.05802929e-03] [-2.58177670e-03 -1.05802929e-03 -9.99999462e-01]] Is it a good fit?2 x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff [[ 7.32686924e-07 4.75306325e-10 -1.70581000e-03] [ 4.75306325e-10 7.32686924e-07 -5.39313733e-04] [-1.70581000e-03 -5.39313733e-04 -9.99999462e-01]] Is it a good fit?2 x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff [[ 5.85509650e-07 4.75306325e-10 -1.36058178e-03] [ 4.75306325e-10 5.85509650e-07 -3.02023162e-04] [-1.36058178e-03 -3.02023162e-04 -9.99999462e-01]] Is it a good fit?2 x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff [[-3.66265272e-07 4.75306325e-10 9.41946218e-04] [ 4.75306325e-10 -3.66265272e-07 6.56238900e-04] [ 9.41946218e-04 6.56238900e-04 -9.99999462e-01]] Is it a good fit?1 x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff [[-4.88584245e-07 4.75306325e-10 1.25433906e-03] [ 4.75306325e-10 -4.88584245e-07 7.83008428e-04] [ 1.25433906e-03 7.83008428e-04 -9.99999462e-01]] Is it a good fit?1 x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff [[-7.02470516e-07 4.75306325e-10 1.80119649e-03] [ 4.75306325e-10 -7.02470516e-07 9.88621422e-04] [ 1.80119649e-03 9.88621422e-04 -9.99999462e-01]] Is it a good fit?1 x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff [[-1.26810745e-06 4.75306325e-10 3.23967294e-03] [ 4.75306325e-10 -1.26810745e-06 1.51545546e-03] [ 3.23967294e-03 1.51545546e-03 -9.99999462e-01]] Is it a good fit?2 x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff [[-3.00324635e-06 4.75306325e-10 7.67149042e-03] [ 4.75306325e-10 -3.00324635e-06 2.95940159e-03] [ 7.67149042e-03 2.95940159e-03 -9.99999462e-01]] Is it a good fit?2 x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff [[ 2.78928085e-06 4.75306325e-10 -7.07337439e-03] [ 4.75306325e-10 2.78928085e-06 -2.04190106e-03] [-7.07337439e-03 -2.04190106e-03 -9.99999462e-01]] Is it a good fit?2 x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff [[ 1.60175355e-06 4.75306325e-10 -4.06730316e-03] [ 4.75306325e-10 1.60175355e-06 -8.32599901e-04] [-4.06730316e-03 -8.32599901e-04 -9.99999462e-01]] Is it a good fit?2 x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff [[-2.66292824e-07 4.75306325e-10 7.36336097e-04] [ 4.75306325e-10 -2.66292824e-07 4.76573016e-04] [ 7.36336097e-04 4.76573016e-04 -9.99999462e-01]] Is it a good fit?1 x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff [[-3.25284852e-07 4.75306325e-10 8.98235317e-04] [ 4.75306325e-10 -3.25284852e-07 5.20749145e-04] [ 8.98235317e-04 5.20749145e-04 -9.99999462e-01]] Is it a good fit?1 x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff [[-4.03746681e-07 4.75306325e-10 1.11518083e-03] [ 4.75306325e-10 -4.03746681e-07 5.67327951e-04] [ 1.11518083e-03 5.67327951e-04 -9.99999462e-01]] Is it a good fit?1 x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff [[-5.21047365e-07 4.75306325e-10 1.43702280e-03] [ 4.75306325e-10 -5.21047365e-07 6.26811405e-04] [ 1.43702280e-03 6.26811405e-04 -9.99999462e-01]] Is it a good fit?1 x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff [[-7.03250637e-07 4.75306325e-10 1.93891868e-03] [ 4.75306325e-10 -7.03250637e-07 6.94528824e-04] [ 1.93891868e-03 6.94528824e-04 -9.99999462e-01]] Is it a good fit?2 x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff [[-9.96262090e-07 4.75306325e-10 2.74601838e-03] [ 4.75306325e-10 -9.96262090e-07 7.62551767e-04] [ 2.74601838e-03 7.62551767e-04 -9.99999462e-01]] Is it a good fit?2 x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff [[-1.36544789e-06 4.75306325e-10 3.75679371e-03] [ 4.75306325e-10 -1.36544789e-06 7.67007660e-04] [ 3.75679371e-03 7.67007660e-04 -9.99999462e-01]] Is it a good fit?2 x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff [[-2.06004470e-07 4.75306325e-10 6.09752267e-04] [ 4.75306325e-10 -2.06004470e-07 3.67305162e-04] [ 6.09752267e-04 3.67305162e-04 -9.99999462e-01]] Is it a good fit?1 x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff [[-2.38580862e-07 4.75306325e-10 7.06047598e-04] [ 4.75306325e-10 -2.38580862e-07 3.80434787e-04] [ 7.06047598e-04 3.80434787e-04 -9.99999462e-01]] Is it a good fit?2 x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff [[-2.77778132e-07 4.75306325e-10 8.22098410e-04] [ 4.75306325e-10 -2.77778132e-07 3.89069642e-04] [ 8.22098410e-04 3.89069642e-04 -9.99999462e-01]] Is it a good fit?1 x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff [[-3.25928641e-07 4.75306325e-10 9.65053368e-04] [ 4.75306325e-10 -3.25928641e-07 3.89726595e-04] [ 9.65053368e-04 3.89726595e-04 -9.99999462e-01]] Is it a good fit?1 x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff [[-3.89253374e-07 4.75306325e-10 1.15176103e-03] [ 4.75306325e-10 -3.89253374e-07 3.81414201e-04] [ 1.15176103e-03 3.81414201e-04 -9.99999462e-01]] Is it a good fit?2 x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff [[-4.61463753e-07 4.75306325e-10 1.36490746e-03] [ 4.75306325e-10 -4.61463753e-07 3.52580456e-04] [ 1.36490746e-03 3.52580456e-04 -9.99999462e-01]] Is it a good fit?2 x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff [[-5.47788559e-07 4.75306325e-10 1.61837109e-03] [ 4.75306325e-10 -5.47788559e-07 2.99621247e-04] [ 1.61837109e-03 2.99621247e-04 -9.99999462e-01]] Is it a good fit?3 x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff [[-1.64699627e-07 4.75306325e-10 5.20141144e-04] [ 4.75306325e-10 -1.64699627e-07 2.94129674e-04] [ 5.20141144e-04 2.94129674e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff [[-1.84993418e-07 4.75306325e-10 5.83879266e-04] [ 4.75306325e-10 -1.84993418e-07 2.95225429e-04] [ 5.83879266e-04 2.95225429e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff [[-2.07921155e-07 4.75306325e-10 6.56266515e-04] [ 4.75306325e-10 -2.07921155e-07 2.90163246e-04] [ 6.56266515e-04 2.90163246e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff [[-2.34890378e-07 4.75306325e-10 7.42205549e-04] [ 4.75306325e-10 -2.34890378e-07 2.78283250e-04] [ 7.42205549e-04 2.78283250e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff [[-2.64584114e-07 4.75306325e-10 8.37927431e-04] [ 4.75306325e-10 -2.64584114e-07 2.55953143e-04] [ 8.37927431e-04 2.55953143e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff [[-3.03670268e-07 4.75306325e-10 9.61501514e-04] [ 4.75306325e-10 -3.03670268e-07 2.22943889e-04] [ 9.61501514e-04 2.22943889e-04 -9.99999462e-01]] Is it a good fit?2 x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff [[-3.31767296e-07 4.75306325e-10 1.05090498e-03] [ 4.75306325e-10 -3.31767296e-07 1.73365250e-04] [ 1.05090498e-03 1.73365250e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff [[-1.45522900e-07 4.75306325e-10 4.26517799e-04] [ 4.75306325e-10 -1.45522900e-07 2.63232521e-04] [ 4.26517799e-04 2.63232521e-04 -9.99999462e-01]] Is it a good fit?5 x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff [[-1.51287709e-07 4.75306325e-10 5.05660579e-04] [ 4.75306325e-10 -1.51287709e-07 2.39860394e-04] [ 5.05660579e-04 2.39860394e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff [[-1.65183726e-07 4.75306325e-10 5.53983800e-04] [ 4.75306325e-10 -1.65183726e-07 2.27513767e-04] [ 5.53983800e-04 2.27513767e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff [[-1.81567497e-07 4.75306325e-10 6.10196583e-04] [ 4.75306325e-10 -1.81567497e-07 2.11818232e-04] [ 6.10196583e-04 2.11818232e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff [[-1.97176154e-07 4.75306325e-10 6.65575727e-04] [ 4.75306325e-10 -1.97176154e-07 1.86834572e-04] [ 6.65575727e-04 1.86834572e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff [[-2.14258821e-07 4.75306325e-10 7.22344515e-04] [ 4.75306325e-10 -2.14258821e-07 1.57804687e-04] [ 7.22344515e-04 1.57804687e-04 -9.99999462e-01]] Is it a good fit?2 x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff [[-1.65240971e-07 4.75306325e-10 5.30559597e-04] [ 4.75306325e-10 -1.65240971e-07 1.66345162e-04] [ 5.30559597e-04 1.66345162e-04 -9.99999462e-01]] Is it a good fit?5
plt.figure(figsize=(15, 10))
center_leds = center_cords_circle[3, 3]
y_unit_len = center_cords_circle[3, 2][1] - center_cords_circle[3, 3][1]
x_unit_len = center_cords_circle[4, 3][0] - center_cords_circle[3, 3][0]
center_leds, y_unit_len, x_unit_len
x_axis = [center_leds[0] - i*x_unit_len for i in range(-4, 4)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-4, 4)]
xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')
x_cords = [x if x != 0.0 else np.nan for x in center_cords_circle.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_circle.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range_circle.reshape(-1)]
plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()
count = 0
for x in range(12, 19):
for y in range(12, 19):
plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
count += 1
plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
axis_lens_pro = axis_lens_circle.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')
xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])
cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')
cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')
plt.show()
axis_lens_pro = axis_lens_circle.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
def fit_ellipse_4_stupid_ransac(img_name, img_8, g_img_8, num_points_to_sample=2000, iterations=10, save_path=None):
""" Find ellipse center position vector where A=-3.24326055e-07, B=0, C=-3.27813554e-07, F=-9.99999462e-01
Parameters
----------
Returns
-------
Ellipse
Ellipse in an image.
"""
ellipsoid_points=[]
img = np.array(g_img_8, dtype=np.uint8)
op = contour_points_v2(img)
op[:10, :] = 0
op[-10:, :] = 0
op[:, :10] = 0
op[:, -10:] = 0
x, y = np.where(op > 0)
final_cnt = np.vstack([y, x]).T
if final_cnt is None:
return None, None, None
else:
if final_cnt.shape[0] > num_points_to_sample:
final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=num_points_to_sample))), axis=0)
# print(final_cnt.shape)
# print(final_cnt)
ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
# print(ellipsoid_points.shape)
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
is_det_ellipse = True
# C = -3.24326055e-07
B = 4.75306325e-10
# A = -3.27813554e-07
F = -9.99999462e-01*2
# b = A*x**2 + B*x*y + C*y**2 + F
x, y = final_cnt[:, 1], final_cnt[:, 0]
b = B*x*y + F
eqs = np.array([x**2 + y**2, x, y]).T
sq_fit = np.linalg.lstsq(eqs, -b, rcond=None)
A, E, D = sq_fit[0]/2
residuals = np.matmul(np.array([x**2 + y**2, x, y]).T, [2*A, 2*E, 2*D]) + b
print(f'0: {np.linalg.norm(residuals)}')
# print(np.sort(residuals)[::-1][:15]**2)
for i in range(iterations):
argmax = residuals.argsort()[-20:][::-1]
final_cnt = np.delete(final_cnt, argmax, axis=0)
x, y = final_cnt[:, 1], final_cnt[:, 0]
b = B*x*y + F
eqs = np.array([x**2 + y**2, x, y]).T
sq_fit = np.linalg.lstsq(eqs, -b, rcond=None)
A, E, D = sq_fit[0]/2
residuals = np.matmul(np.array([x**2 + y**2, x, y]).T, [2*A, 2*E, 2*D]) + b
print(f'{i+1}: {np.linalg.norm(residuals)}')
# print(np.sort(residuals)[::-1][:15]**2)
matC = np.array([[A, B, D/2], [B, A, E/2], [D/2, E/2, F/2]])
E = Ellipse( matC )
# print(Ellipse.is_ellipsoid( E.C ))
# print(matC)
# print(np.linalg.lstsq(eqs, -b, rcond=None))
if Ellipse.is_ellipsoid( E.C ):
S,R,T = Ellipse.decompose_ellipsoid(E.C)
tx = T[0,2]
ty = T[1,2]
smin = S[0,0]
smaj = S[1,1]
theta = np.arctan2(R[1,0], R[0, 0])
# print(theta)
# f = E.C[-1, -1]
# print(E.C)
else:
is_det_ellipse = False
if is_det_ellipse:
for i in range( 0, ellipsoid_points.shape[1] ):
draw_cross(img_8, (ellipsoid_points[0,i],ellipsoid_points[1,i]), (0, 0, 255), length=20, thickness=5)
cv2.ellipse(img_8, (int(tx), int(ty)), (int(smin), int(smaj)), theta, 0, 360, (255, 0, 0), 5)
draw_cross(img_8, (tx, ty), (0, 255, 0))
# plt.imshow(img_8)
# plt.show()
cv2.imshow('Resized Image', cv2.resize(img_8, (1080, 720), interpolation = cv2.INTER_LINEAR))
cv2.waitKey(0)
cv2.destroyAllWindows()
if save_path is not None:
cv2.imwrite(f'{save_path}/{img_name[:-5]}.png', img_8)
return (int(tx), int(ty)), (int(smin), int(smaj)), theta
else:
return None, None, None
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff"
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
# np.random.seed(2)
# # grid based on thumbnail images
# fit_range_circle_iter = np.zeros([7, 7], dtype=np.uint8)
# center_cords_circle_iter = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_circle_iter = np.zeros([7, 7, 2], dtype=np.float64)
# save_dir = "C:/Users/Kazim/Desktop/Overall/VignettingCircleDetection/310823_circle_iter_fit_tiff_size"
# for x in range(12, 19):
# for y in range(12, 19):
# dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
# file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
# print('x:', x, 'y:', y, file_path)
# img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
# _, g_img, _ = cv2.split(img)
# img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
# g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
# center, axis_size, f = fit_ellipse_4_stupid_ransac(file_path, img_8, g_img_8, num_points_to_sample=2000,
# iterations=10, save_path=save_dir)
# if center is not None:
# good_fit = input("Is it a good fit?")
# fit_range_circle_iter[x-12, y-12] = int(good_fit)
# center_cords_circle_iter[x-12, y-12] = center
# axis_lens_circle_iter[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff 0: 20.5704220924657 1: 20.411638772482917 2: 20.143170477624405 3: 19.643555622764538 4: 18.157343271262192 5: 12.29152406008125 6: 1.6314344090668467 7: 1.5779969524814406 8: 1.5397842754725057 9: 1.5057562645965468 10: 1.4750242120243569 Is it a good fit?1 x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff 0: 5.8653172341717195 1: 5.582038242007653 2: 5.290316054749441 3: 4.99524088877848 4: 4.688428229662138 5: 4.384281332500998 6: 4.159156402148534 7: 3.9919146573340787 8: 3.896977109982126 9: 3.8100337626772114 10: 3.728288447442476 Is it a good fit?2 x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff 0: 3.5902079940248495 1: 3.4591804994062274 2: 3.352404816533077 3: 3.252449759584019 4: 3.1566990197310405 5: 3.0652135186672074 6: 2.9752301506202152 7: 2.88882638793891 8: 2.805712317154255 9: 2.726922709902689 10: 2.6580234687286586 Is it a good fit?1 x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff 0: 1.2546587604834094 1: 1.1542997949086367 2: 1.1294327033060079 3: 1.1066635156415723 4: 1.084783180998509 5: 1.0633381977633445 6: 1.042573035707542 7: 1.0220822389276405 8: 1.0017416461388435 9: 0.9817834636235067 10: 0.9619089317818343 Is it a good fit?2 x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff 0: 1.6018928596932482 1: 1.5664494426892628 2: 1.5339115849704665 3: 1.500744497892927 4: 1.4653539103343924 5: 1.428099202573096 6: 1.3887127843657292 7: 1.3462691099820296 8: 1.3001138069724572 9: 1.2499513509742812 10: 1.1942545077022388 Is it a good fit?2 x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff 0: 0.9840795107478418 1: 0.9712018450403989 2: 0.9589012605095123 3: 0.9470017227369132 4: 0.9348025946912554 5: 0.9222568229441381 6: 0.9093369119286429 7: 0.8957977916671148 8: 0.8815130978766169 9: 0.8660658817610803 10: 0.8492869288375672 Is it a good fit?2 x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff 0: 5.313513662083277 1: 5.2958757624450765 2: 5.284872156015689 3: 5.2774754825143955 4: 5.273167276780472 5: 5.270584038926708 6: 5.268300161837872 7: 5.26608226857869 8: 5.263856011537344 9: 5.261604854370673 10: 5.259313119511193 Is it a good fit?2 x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff 0: 0.3348373587586304 1: 0.3198855210954769 2: 0.30857885200760826 3: 0.2992491989086529 4: 0.2907914218059756 5: 0.2829356881632068 6: 0.2759106175425677 7: 0.26939994331511985 8: 0.263404992853191 9: 0.2577527599914932 10: 0.2523856228249853 Is it a good fit?2 x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff 0: 1.133830369580832 1: 1.0913247217060837 2: 1.0595968195145298 3: 1.0364629845154034 4: 1.0166541339175221 5: 0.99992446490275 6: 0.9843870199978809 7: 0.9690844687558989 8: 0.9543601446777084 9: 0.9399318415421593 10: 0.925669305251249 Is it a good fit?1 x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff 0: 10.014975678733723 1: 8.736813312473847 2: 7.285188742355118 3: 5.585806810610128 4: 3.3550069448170623 5: 2.326445067428788 6: 2.2450556671417554 7: 2.17560762056779 8: 2.1143668200869246 9: 2.0577821972815817 10: 2.005698352312667 Is it a good fit?2 x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff 0: 8.449862487866428 1: 8.281831187514262 2: 8.111563298705253 3: 7.933473262248377 4: 7.74480478591166 5: 7.5440945167361315 6: 7.325519591469758 7: 7.085659152459896 8: 6.827696309282607 9: 6.62957859472071 10: 6.4595303995783615 Is it a good fit?2 x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff 0: 3.5415542098585013 1: 3.502860764161347 2: 3.4641045461603373 3: 3.4244660793408563 4: 3.383189937510767 5: 3.3404744132414 6: 3.29611459948416 7: 3.2500554535400417 8: 3.201492238725641 9: 3.1498985314513255 10: 3.0951824092089657 Is it a good fit?2 x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff 0: 1.6038564087150302 1: 1.583109262360901 2: 1.5641762131799806 3: 1.545549825180772 4: 1.527333805872587 5: 1.5098587662065224 6: 1.4931290081007116 7: 1.477108409680803 8: 1.4615903350510406 9: 1.4461928972308806 10: 1.4308435474583288 Is it a good fit?2 x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff 0: 1.1614312260064732 1: 1.1403060811506096 2: 1.1218188209451867 3: 1.1045587254619482 4: 1.0887381293694771 5: 1.0742570136851255 6: 1.0612811566690914 7: 1.0495538179848345 8: 1.0383953016451244 9: 1.027573108274712 10: 1.01736698286472 Is it a good fit?3 x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff 0: 0.2377662069226674 1: 0.2007037735688992 2: 0.1824014029687429 3: 0.17501977185692028 4: 0.16956335068036735 5: 0.16468185005782088 6: 0.16019835303040802 7: 0.15600696054055133 8: 0.15201260659908883 9: 0.14813088017330425 10: 0.14443523170867276 Is it a good fit?1 x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff 0: 0.3577055618340971 1: 0.3219689865115248 2: 0.29064384223837686 3: 0.265887847683449 4: 0.2516295875334015 5: 0.2420086532067681 6: 0.23346940044576006 7: 0.22555513530921353 8: 0.21797578686977662 9: 0.2107410898751867 10: 0.20414614341215606 Is it a good fit?1 x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff 0: 1.0550361599683238 1: 0.9991159235416885 2: 0.9470189788145269 3: 0.8971464233064096 4: 0.8525983397011151 5: 0.8177678097397628 6: 0.7865985295488093 7: 0.7578168877488449 8: 0.7318163495359703 9: 0.7070716887966484 10: 0.6835805454215517 Is it a good fit?1 x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff 0: 2.4568239174465947 1: 2.2719677283256763 2: 2.092149348106298 3: 1.926709042519453 4: 1.8289223537232704 5: 1.7374786256789276 6: 1.647378178585114 7: 1.5610120401203935 8: 1.477782897717219 9: 1.3983857792268413 10: 1.3286610079822045 Is it a good fit?2 x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff 0: 10.020378807838993 1: 9.666403292162503 2: 9.370770979185151 3: 9.118084436087281 4: 8.887694680076471 5: 8.676626072158925 6: 8.471285058065924 7: 8.262413970799477 8: 8.04740387022261 9: 7.828037707476127 10: 7.602018984370476 Is it a good fit?2 x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff 0: 11.635329947331929 1: 11.564397336490309 2: 11.495785899396534 3: 11.427598982545883 4: 11.359095513258142 5: 11.289593664813564 6: 11.219139646748946 7: 11.146584483309633 8: 11.073178878804866 9: 10.999123788165253 10: 10.924780484242858 Is it a good fit?3 x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff 0: 5.616936156093774 1: 5.574257972431561 2: 5.542292585710326 3: 5.512297678750731 4: 5.482770235093951 5: 5.453972712411235 6: 5.42535145086847 7: 5.3973489914734865 8: 5.36929856498738 9: 5.341149018662106 10: 5.312493774666667 Is it a good fit?3 x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff 0: 0.3013473542410165 1: 0.29053700798804216 2: 0.28304285660288603 3: 0.2764010788084259 4: 0.26992084529138294 5: 0.26347240296586566 6: 0.2568317448736444 7: 0.25032315128171456 8: 0.24375464931444205 9: 0.23672457143180273 10: 0.22901635734512021 Is it a good fit?1 x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff 0: 0.42492208666900233 1: 0.41191610924434496 2: 0.40066777012826765 3: 0.39057113804152765 4: 0.38116935048900646 5: 0.3720623900701497 6: 0.36349998780586773 7: 0.3553431415001986 8: 0.3476405535728397 9: 0.34067371723416795 10: 0.333991772738715 Is it a good fit?1 x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff 0: 0.39989671490546497 1: 0.3783492190819314 2: 0.36147574741447874 3: 0.3483905653557774 4: 0.33799711650300535 5: 0.32873046772251185 6: 0.3203720219284474 7: 0.31248841718301756 8: 0.3049908648747955 9: 0.2977709594563425 10: 0.29080241544475777 Is it a good fit?1 x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff 0: 2.0292570771141105 1: 1.8191418502815577 2: 1.6039281207342626 3: 1.3847949550976428 4: 1.1774467813600529 5: 0.9799803338126872 6: 0.9190078058439943 7: 0.8716444956894469 8: 0.8283253702370312 9: 0.7900285219489808 10: 0.7580279437566125 Is it a good fit?1 x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff 0: 2.372760657014323 1: 2.309355794813653 2: 2.2557655129457825 3: 2.2073217712323303 4: 2.161690904703809 5: 2.116305573168592 6: 2.069990715554002 7: 2.02351681011232 8: 1.9761004745861885 9: 1.9286433096262847 10: 1.8808164720321137 Is it a good fit?2 x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff 0: 4.422923530579537 1: 4.301079526438787 2: 4.1931284997878935 3: 4.0629489890794375 4: 3.8975941057304166 5: 3.6956071775516937 6: 3.4669724488683986 7: 3.281476693133121 8: 3.1331534541086836 9: 2.94223186108211 10: 2.7300502542821565 Is it a good fit?4 x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff 0: 7.2196186072880915 1: 7.089296329189396 2: 6.961239330430815 3: 6.832470992590163 4: 6.702381550149271 5: 6.572596090852814 6: 6.44158949193341 7: 6.314016416593842 8: 6.201068447999578 9: 6.103637758183448 10: 6.007688289221982 Is it a good fit?3 x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff 0: 0.31224661233731027 1: 0.3026496129800401 2: 0.2954287031732409 3: 0.28866427829588126 4: 0.28162538274567234 5: 0.27412770149942467 6: 0.2662211792487556 7: 0.2576046267231648 8: 0.24783530221649686 9: 0.23667148596132373 10: 0.22363205394300517 Is it a good fit?2 x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff 0: 0.2813841061870855 1: 0.2747195004822052 2: 0.2690200966981868 3: 0.2636957863235934 4: 0.2586181077977398 5: 0.2536405894039849 6: 0.24869718022169954 7: 0.24380390643065533 8: 0.23879822387679253 9: 0.23360215975033577 10: 0.22814620944988975 Is it a good fit?1 x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff 0: 0.231990751344788 1: 0.22448018515277954 2: 0.2180389042944469 3: 0.2122511915927251 4: 0.20722233594279652 5: 0.20262011059956922 6: 0.19825858917639597 7: 0.19402234225394563 8: 0.18988499633295872 9: 0.18584610478200797 10: 0.1818792597833647 Is it a good fit?1 x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff 0: 0.7188694709796103 1: 0.6606489243642706 2: 0.6148027142262953 3: 0.5750277436602184 4: 0.5381782847063585 5: 0.5018953757971297 6: 0.46574984979784445 7: 0.4327739652028631 8: 0.403608118736117 9: 0.3838120733474389 10: 0.3701203897753593 Is it a good fit?1 x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff 0: 1.4089737833576692 1: 1.3788941478343097 2: 1.3508803538022929 3: 1.3241038093417132 4: 1.2995112748619726 5: 1.2764975144428465 6: 1.2542006737279232 7: 1.2313468095344453 8: 1.2081187361522066 9: 1.1846076137990784 10: 1.1604339988780965 Is it a good fit?2 x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff 0: 2.1051620110250493 1: 2.0262136777185766 2: 1.9621093044493703 3: 1.9029866354180953 4: 1.8492383165298802 5: 1.78555414579459 6: 1.7044910993733706 7: 1.602139665019008 8: 1.4764937925620611 9: 1.2980751165178375 10: 1.0847245901700924 Is it a good fit?3 x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff 0: 2.6255241091914647 1: 2.4987195132007485 2: 2.420327983806968 3: 2.332389465849195 4: 2.236634311232568 5: 2.1348584985469037 6: 2.019398743482966 7: 1.872816648507637 8: 1.6908061244024348 9: 1.4803080707012042 10: 1.3274829039897278 Is it a good fit?4 x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff 0: 0.13401237965073892 1: 0.129432596935074 2: 0.1252631602301332 3: 0.12166853638288813 4: 0.11846101571754876 5: 0.11546745898400519 6: 0.11283047092420898 7: 0.11039949415248872 8: 0.10811695501844279 9: 0.10590602972671448 10: 0.10383135241649526 Is it a good fit?2 x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff 0: 0.12738582820914102 1: 0.11744208829963568 2: 0.11117452044799904 3: 0.10618860593791832 4: 0.10204328360799371 5: 0.0986132134697864 6: 0.09542294697924678 7: 0.09249793468460811 8: 0.08978666433709857 9: 0.08728001378582963 10: 0.08492272313688032 Is it a good fit?2 x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff 0: 0.19347801846232482 1: 0.18379865795646033 2: 0.17535667802098656 3: 0.16888238697261584 4: 0.16349953820766422 5: 0.15870617959796388 6: 0.15420743673511042 7: 0.1499337436482179 8: 0.14573921286238053 9: 0.14183025197627988 10: 0.1382034892821991 Is it a good fit?1 x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff 0: 0.6544198933025352 1: 0.6331518293452286 2: 0.6154835776599894 3: 0.6073227010628626 4: 0.5995118819158949 5: 0.5915891649417299 6: 0.5835240296814301 7: 0.5752201273230771 8: 0.5667649106762465 9: 0.5581074304268094 10: 0.5492288632333817 Is it a good fit?2 x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff 0: 0.8534051290794373 1: 0.828214919595704 2: 0.8005796405504847 3: 0.770351259979737 4: 0.7468177679544643 5: 0.727263441650287 6: 0.7038323768768756 7: 0.6735652994382467 8: 0.633398692442303 9: 0.5772576867633307 10: 0.5126497538653242 Is it a good fit?3 x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff 0: 0.9114550493011583 1: 0.7916623322151255 2: 0.648339110569025 3: 0.5330626310585324 4: 0.48415823410722947 5: 0.46007797409208456 6: 0.43500750864175863 7: 0.41250208800183313 8: 0.39599552535819105 9: 0.3853444519565576 10: 0.3760450998105657 Is it a good fit?3 x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff 0: 1.0199021295951793 1: 0.9297353328997268 2: 0.8606242055516888 3: 0.7877509888542633 4: 0.7067557733718133 5: 0.6157466159337168 6: 0.525440611537567 7: 0.4705458155497611 8: 0.4424962559401255 9: 0.4160009042775519 10: 0.3899266002836627 Is it a good fit?2 x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff 0: 12.372346792013559 1: 11.39455880780588 2: 10.226206395761132 3: 8.946495839851549 4: 7.222329463849575 5: 6.015753874149101 6: 4.272571679770303 7: 3.055794606125172 8: 1.9407424440836454 9: 0.1254506683049625 10: 0.1191479074715457 Is it a good fit?2 x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff 0: 0.15051619338828395 1: 0.14542505203572068 2: 0.1411964048866287 3: 0.13745409892633118 4: 0.13405979461103962 5: 0.1307883180364522 6: 0.12711917606080092 7: 0.12055244718917445 8: 0.11769407808175629 9: 0.11498882108355644 10: 0.11242827100649683 Is it a good fit?3 x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff 0: 0.24453599651961 1: 0.2391690865772408 2: 0.23411270874244183 3: 0.22916944878026174 4: 0.22443448857057716 5: 0.21965734447928972 6: 0.21497105329403718 7: 0.2106978589025084 8: 0.20671815945840336 9: 0.20291825353474194 10: 0.19910562748428617 Is it a good fit?2 x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff 0: 0.4797979493962809 1: 0.4659191891821974 2: 0.45368511272072615 3: 0.44073771079022594 4: 0.42649540274460174 5: 0.40980421125264344 6: 0.39025625449679846 7: 0.3749944252341886 8: 0.3654218867146383 9: 0.35597015662365883 10: 0.34640391378428476 Is it a good fit?2 x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff 0: 0.6019285018047863 1: 0.5650090332598229 2: 0.5218588500302644 3: 0.48222051420463025 4: 0.46023697736640273 5: 0.4400021711815074 6: 0.41608810528694395 7: 0.386206859723485 8: 0.3521841161300144 9: 0.32053888114460133 10: 0.2859715074166922 Is it a good fit?3 x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff 0: 0.6066184611995026 1: 0.5771131622850886 2: 0.5450949549683325 3: 0.5103371745038922 4: 0.4836031240352168 5: 0.45783203543063866 6: 0.43308624175285815 7: 0.4098906352571529 8: 0.3868128561751656 9: 0.36454392634747496 10: 0.3416481105895158 Is it a good fit?2 x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff 0: 14.142003659512072 1: 12.626426154055927 2: 9.616114276770336 3: 0.616739113038089 4: 0.5838648610653595 5: 0.5557056913369234 6: 0.5286274090762597 7: 0.5027140147309612 8: 0.47672588443512115 9: 0.45126333232542465 10: 0.4262611148576096 Is it a good fit?2
plt.figure(figsize=(15, 10))
center_leds = center_cords_circle_iter[3, 3]
y_unit_len = center_cords_circle_iter[3, 2][1] - center_cords_circle_iter[3, 3][1]
x_unit_len = center_cords_circle_iter[4, 3][0] - center_cords_circle_iter[3, 3][0]
center_leds, y_unit_len, x_unit_len
x_axis = [center_leds[0] - i*x_unit_len for i in range(-4, 4)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-4, 4)]
xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')
x_cords = [x if x != 0.0 else np.nan for x in center_cords_circle_iter.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_circle_iter.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range_circle_iter.reshape(-1)]
plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()
count = 0
for x in range(12, 19):
for y in range(12, 19):
plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
count += 1
plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
axis_lens_pro = axis_lens_circle_iter.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')
xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])
cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')
cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none',
extent=[12, 19, 12, 19]
)
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')
plt.show()
axis_lens_pro = axis_lens_circle_iter.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan
plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
